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By computer simulations of systems of ellipsoids, we study the influence of the isotropic/nematic 
phase transition on the direct correlation functions (DCF) in anisotropic fluids. The DCF is de- 
termined from the pair distribution function by solving the full Ornstein-Zernike equation, without 
any approximations. Using a suitable molecular-fixed reference frame, we can distinguish between 
two qualitatively different contributions to the DCF: One which preserves rotational invariance, and 
one which breaks it and vanishes in the isotropic phase. We find that the symmetry preserving 
contribution is barely affected by the phase transition. However, symmetry breaking contributions 
emerge in the nematic phase and may become quite substantial. Thus the DCF in a nematic fluid is 
not rotationally invariant. In the isotropic fluid, the DCF is in good agreement with the prediction 
of the Percus-Yevick theory. 



I. INTRODUCTION 
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The pair direct correlation function (DCF) is a central 
quantity in the theory of liquids. Through the Ornstein- 
Zernike relation 1 , it is related to the pair distribution 
function, which is accessible experimentally. Compared 
to the latter, it often has a much simpler structure, be- 
cause the Ornstein-Zernike relation eliminates to a large 
extent the contributions of third particles to the pair cor- 
relations. Moreover, the DCF plays a key role in den- 
sity functional theories, since it is the second functional 
derivative of the excess free energy with respect to the 
local density 1,2 . Much effort has therefore been devoted 
to the investigation of DCFs in molecular fluids. 

Most studies have considered isotropic fluids: Early in- 
tegral equation approaches are due to Wulf 7 , and Chen 
and Steele 6 . Chen and Steele generalized the Ornstein- 
Zernike equation to the case of particles with orientation 
dependent interactions 6,31 and calculated the structure of 
diatomic fluids using a Percus-Yevick 4 closure. Percus- 
Yevick and Hypernetted Chain 1 theories for isotropic 
fluids with orientation dependent interactions have also 
been introduced by Pynn 5 , Fries and Patek 8 , and Singh 
and coworkers 9 . Several authors have calculated DCFs in 
isotropic liquid crystals using theories of this kind 5 ' 10-13 . 
The resulting pair correlation functions were generally in 
reasonable agreement with available simulation data. 

Alternative approaches have been pursued, e. g., by 
Wulf 7,14 , Rickayzen and coworkers 15-17 , and Chamoux 
and Perera 18 . Rickayzen et al choose a geometrically 
motivated Ansatz for the form of the DCF and optimize 
the parameters such that the Percus-Yevick relation is 
fulfilled in a given limit. Chamoux and Perera derive 
approximations for the DCF from the Rosenfcld den- 
sity functional 19 . For fluids of elongated particles, their 
DCF is similar to that obtained with the Percus-Yevick 
or Hypernetted chain calculations of Perera et al 10 . A 
number of other authors have suggested simple approxi- 
mate expressions for the DCF in isotropic molecular flu- 
ids 5 ' 14 ' 21-24 . Some of these fit simulation data surpris- 



ingly well 25 . The direct determination of the DCF from 
computer simulations is rather tedious. Allen et al 25 have 
presented a method which allows to calculate the DCF 
from pair correlation data, and applied it to study the 
DCF in fluids of ellipsoids 25 
cles 27,28 . 
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The DCF has thus been studied quite intensely in 
isotropic fluids. In contrast, only few studies have con- 
sidered the DCF in anisotropic, e. g., nematic fluids. In 
a nematic liquid crystal, the particles remain position- 
ally disordered, but align preferentially along one (ar- 
bitrary) direction 29,30 . The isotropy of space is spon- 
taneously broken. Consequently, the pair correlation 
functions lose their rotational invariance. Moreover, soft 
Goldstone modes must be present, since the orientational 
order breaks a continuous symmetry: The structure fac- 
tor diverges in the k — > limit, and the pair distribution 
function exhibits a long range elastic tail. This implies 
that the DCF must fulfill certain conditions in the k — > 
limit, which have been derived by Gubbins 31 , and later 
in a different manner by Zhong and Petschek 32 . 

Workman and Fixman 33 have formulated the gen- 
eral form of the Ornstein-Zernike equation for the case 
of anisotropic fluids. They also suggested a Percus- 
Yevick closure which is based on a density functional 
about an isotropic reference state, but can be applied 
to both isotropic and anisotropic liquids. In fact, the 
Percus-Yevick and the Hypernetted chain closure lend 
themselves to a rather straightforward generalization for 
anisotropic fluids 34 . Caillol and coworkers 34 have used 
these closures to calculate the structure of perfectly 
aligned fluids. Zhong and Petschek 32 have analyzed the 
diagrammatic expansion of the DCF in this approach, 
and proved that at least the Percus-Yevick closure must 
fail to reproduce the soft Goldstone modes in the gen- 
eral case of spontaneous partial order. In order to fix the 
problem without resorting to a reference state, they pro- 
posed a modified Percus-Yevick closure which premises 
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that the DCF is rotationally invariant. Given that the 
DCF is vaguely related to an "effective interaction poten- 
tial" between particles 4 ' 32 , this assumption seems plau- 
sible. Nevertheless, a clearcut test of the hypothesis is 
clearly desirable. 

In practice, theoretical studies of DCFs in partially 
ordered nematic fluids have mostly been concerned with 
systems with separable interactions 35 ' 36 , where the inter- 
particle potential depends on the spatial separation and 
the orientation of particles independently. Very little is 
still known on the DCF in general anisotropic fluids. Sim- 
ulation studies have been performed by Stelzer et al 37 
(nematic Gay-Berne fluids), and by Zakharov and Ma- 
liniak 38 (5CB molecules). These authors used an "unori- 
ented nematic approximation", which replaces the pair 
correlations by their rotational averages. The simplifi- 
cation has later been questioned by Longa et al 39 . In- 
deed, the elastic properties of the nematic fluid, which 
can be calculated from the DCF using the Poniewierski- 
Stecki equations 40 , don't seem to be captured very well 41 . 
Moreover, an analysis which uses such an approximation 
is clearly not suited to elicit whether or not the DCF 
depends on the direction of prefered alignment in the ne- 
matic fluid. 

In a previous paper 42 , we have presented a method for 
determining the DCF in uniaxial nematic fluids without 
approximations from computer simulations. As a test of 
the method, we have applied it to a nematic fluid of soft 
ellipsoids and calculated the elastic constants from the 
DCF via the Poniewierski-Stecki equations 40 . The re- 
sults were in good agreement with those obtained inde- 
pendently from an analysis of order tensor fluctuations, 
following a procedure by Allen et al 41 . Thus we have 
shown that the method takes into account the elastic 
properties of the fluid in an adequate way. 

In the present work, we use our method to analyze the 
form of the DCF in detail for different state points in the 
nematic and the isotropic phase. In particular, we dis- 
cuss the symmetry properties and the range of the DCF, 
the pair distribution function, and the total correlation 
function. We show that the DCF is truly short ranged 
as expected. The elastic tails and the near-critical long 
range fluctuations close to the nematic-isotropic transi- 
tion disappear. On the other hand, our data show that 
the DCF does reflect the broken symmetry of the nematic 
phase, and the assumption of rotational invariance is not 
correct in our system. 

The paper is organized as follows: In the next section, 
we define the simulation model, introduce the pair cor- 
relation functions and explain how we perform the data 
analysis. The results are presented and discussed in sec- 
tion III. In the isotropic phase, we compare them with 
the prediction of the Percus-Yevick theory. Some com- 
ments on the Percus-Yevick solution and how it breaks 
down at large densities are added in the appendix. We 
summarize and conclude in section IV. 



II. MODEL AND METHOD 

A. simulation model 

We have studied a fluid of soft ellipsoidal particles with 
repulsive pair interactions 
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Xf 2 > 1/2 
otherwise 



(1) 



The function X i2 = <Jo/( r i2 — &12 + &o) depends on the 
distance n 2 between particles 1 and 2, and on the shape 
function 43 



(Ui ■ fi2 + U 2 ■ fl 2 ) 2 
1 + YUi • U 2 

-1/2 

(2) 



(712 (Ui, U 2 ,f 12 ) = cr | 4 ~f 

(ui ■ f 12 - u 2 ■ rig) 

1 - YUi • U 2 



which approximates the contact distance between two el- 
lipsoids of elongation n = \J {I + x)/(l — X) with orien- 
tations Ui and u 2 and center-center vector ri 2 pointing in 
the direction f 12 = r 12 /r 12 . The parameters were k = 3 
and ksT = 0.5eo- 

Four state points were considered, two at densities 
g = 0.24/(7q and 0.283/ctq in the isotropic phase, and 
two at densities 0.287/co an d Q = 0.3/(7o m the ne- 
matic phase. The system size was N = 1000 particles 
at g = 0.24/ag, N = 4000 at g = 0.283/ct 3 , and 0.287/crg, 
and N = 1000 - 8000 at g = 0.3/er 3 ,. The systems with 
the three lower densities were studied in the canonical 
ensemble by Monte Carlo simulation. The data for the 
larger systems (N > 4000) at the density p = 0.3a -1 
were originally produced by G. Germano 42 . They were 
simulated in the microcanonical ensemble using a parallel 
molecular dynamics program on a CRAY T3E. Here the 
moment of inertia of a particle was chosen / = 2.5mo(7o 2 
and the time step At = 0.003<ro y/rnoj^o, where mo is the 
mass of one particle. The results did not depend on the 
simulation method. Run length were 5-10 million MC or 
MD steps, depending on the system size, and data were 
collected every 1000 or 10000 steps. 

Throughout this paper, the units will be defined in 
terms of the energy unit eo, the mass unit mo, the tem- 
perature unit eo/fce, and the length unit ctq. 

The orientational order is characterized as usual by the 
order tensor 29 



JV 



i=l 



(3) 



where the sum i runs over all N particles of orientation 
I is the unit matrix, <g) the dyadic vector product, (•) 
denotes thermal averages. The largest eigenvalue of Q is 
the nematic order parameter VP2, and the correspond- 
ing eigenvector n is the director, which points along the 
direction of prefered alignment. In our case, we have 
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P 2 = 0.47 at g = 0.287/crg, and P 2 = 0.69 at g = 0.3/er 3 ,. 
The densities g = 0.283/og and g = 0.287/og are close 
to the coexisting densities at the isotropic/nematic phase 
transition. Previous constant pressure simulations of the 
same system 44 showed that the transition is only very 
weakly first order, and the width of the coexistence gap 
is of the order of 0.05/<7g. 



B. Data analysis 

Before describing our way to analyze the data, we re- 
call some basic definitions: In a system of particles i with 
orientation and position r^, the one-particle distribu- 
tion can be written as 1 



p( 1 )(u,r) = (^^(u-u i )«5(r-r i )). 



(4) 



and the pair distribution as 

p (2) (ui,ri,u 2 ,r 2 ) (5) 
= (y] Sjm - u i )S(r 1 - Ti)6(u 2 - Uj)5(r 2 - rj)). 

In homogeneous nematic fluids, the pair distribution de- 
pends on the difference vector r\ 2 = ri— r 2 rather than on 
the individual positions ri and r 2 , and the single-particle 
distribution depends only on u (in fact, on |u • n|). In 
that case, the total correlation function /i(m, u 2 , ri 2 ) is 
defined by 



h(ui,U2,ri 2 ) 



/P (2) (llliU2,ri2) 

P W(m)pW(u 2 ) 



1, 



(6) 



and the direct correlation function c(ui, u 2 , ri 2 ) is deter- 
mined by the Ornstein-Zernike equation 1,3 ' 33 

/i(ui,u 2 ,r 12 ) = c(ui,u 2 ,r 12 )+ 

c(ui,u 3 ,ri 3 ) /5 (1) (u 3 ) /i(u3,u 2 ,r 32 )du 3 dr3. (7) 



We choose a coordinate frame where the z axis points 
in the direction of the director n (director frame) 45 . All 
orientation dependent functions are expanded in spheri- 
cal harmonics YJ m (u). This yields 



(8) 



with the bulk number density g, and 



F(ui,u 2 ,r)= ^2 F hmi i 2m2lm (r) (9) 

^mi(ui) ^ 2 m 2 (u 2 ) Yim(f), 



!l,! 2 ,! 



where F stands for p^ 2 \ h or c, and f denotes the unit 
vector r/r. In uniaxially symmetric phases, only real co- 
efficients with mi + TO 2 + m = and even l\ + Z 2 + Z 



enter the expansion. Since our particles have uniaxial 
symmetry, every single Zj is even as well. 

The expansion coefficients of the pair distribution func- 
tion can be determined directly from the simulations 46 



(Y l * imi (u 1 )Y l * 2m2 (u 2 )Y l * m (v))sr, 



(10) 



where the average (-)sr is performed over all molecule 
pairs at distances |ri — r 2 | G [r, r + Sr], and g(r) is the ra- 
dial distribution function, i. e., the average total number 
of molecule pairs divided by N4Trgr 2 5r. The coordinate 
frame was determined separately in each configuration, 
such that the z axis is given by the direction of the di- 
rector (director frame). In general, we have determined 
coefficients for values of Z, k up to Z max = 6. For the high- 
est density 42 , we checked in small systems that the trun- 
cation is sufficient by comparing the results with those 
obtained using Z max = 8. The bin size was Sr = 0.04cr 
and only ranges up to r max = 0.4L [L being the size of 
the simulation box) were considered in order to reduce 
boundary effects 47 . 

The total correlation function h is obtained from the 
spherical harmonics version of eqn. (6), which is simply 
a matrix equation: 

PhLll2m 2 l,n( r ) = ^ ( V^f h f h 8 mi0 6 m2 oS W S m0 (H) 

+ ^l[m 1 l' 2 m 2 lm( r )fl'^ fl' 2 ^ mi mla^m 2 mla ) ' 

ii ui ,1 in 

'l'l ''2> l 2 

with 

Iw m « = J duY^{u)Y Vtmt {u)Y lH>m .,{u) (12) 

The DCF c is most easily calculated in Fourier rep- 
resentation. The expansion coefficients of a function 
F(ui,u 2 ,r) in real space are related to their counter- 
parts in Fourier space by the Hankel transform 3 

i>OC 

Fi imi i 2 rn 2 im(k) = 4ttz' / r 2 ji(kr) F hmi i 2 . m2 i m (r) dr (13) 
Jo 



l\m\l 2 m 2 lm 



(r) 



4tt(- 



(2^) 3 Jo 



k ji(kr)Fi {k)dk (14) 



where ji(kr) is the spherical Bessel function. The 
Ornstein-Zernike equation (7) in Fourier space and spher- 
ical harmonics representation reads 

^l\ra\l 2 m 2 lm 

(k) (k) 

~^~Q ^ ' QimiismsJ'm' (&) hl' 3 m 3 l 2 m 2 l" m" (^) 



•y f,ii(—~\\ m3 V 11 ' 1 " Y l s l 3 l 3 

^Jl 3 \ L ) L mm'm" L m 3 m 3 0' \ 10 / 

which is again a matrix equation. 
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Special care must be taken when performing the Han- 
kel transform (13) in the nematic state: Due to the elas- 
ticity of the nematic fluid, certain coefficients of h (in 
particular those with m l7 m 2 = ±1) have pronounced 
long-range tails and exhibit a 1/r behavior at large dis- 
tances r. Furthermore, large-wavelength fluctuations are 
suppressed in finite systems, and the director frame fol- 
lows the local order if evaluated in a small box. As a 
result, some coefficients (those with I = 0, rrij = ±1) are 
shifted towards zero in finite systems by an amount which 
is almost constant beyond r = o-q (see Figure 1). Before 
performing the Hankcl transform, we thus fit the data 
points beyond distances r > tq to a power law a + b/r, 
shift h(r) by a, and extrapolate it to infinity. At the 
highest density g = 0.3/cTq, we have simulated systems 
of different size (N = 1000, N = 4000, and N = 8000) 
and checked that finite size effects are basically elimi- 
nated by our procedure. The parameter ro was chosen 
r = 4.0/cto at TV = 4000, r = 2.8/a at TV = 1000, and 
r = 5.3/ctq at N = 8000. 
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FIG. 1. Expansion coefficient /1212-100 of the total correla- 
tion function vs. 1/r for systems of different size at density 
p = 0.3o"q~ 3 . Before performing a Hankel transform of this 
curve, the tails are fitted with a function a + b/r, and the 
whole curve is shifted by a. Inset shows the shifted data. The 
curves for the two larger systems are then almost on top of 
each other. 



The data manipulations which finally yield the DCF 
are altogether quite extensive. An independent test of the 
results which checks the procedure is clearly desirable. 
For example, one can calculate the elastic constants from 
the DCF using the Poniewierski-Stecki relations 40 , and 
compare them with values obtained independently from 
the structure factor. As we have reported in our earlier 
paper 42 , the agreement is quite satisfactory. As another 
test, we have verified the validity of a relation originally 
derived by Gubbins 31 , 



p (1)/ K«) 



Ul,a = J C(ui, U 2 , k)| fc=0 P W '{u 2 , z )u 2 , a d\l 2 . 

(16) 



In Fourier space and spherical harmonics representation, 
this reads 



y/W+7) / ln(p( 1 )(u))y ;o (u) du 

= yyi'P + l)/l'Qij'-ioo(fc)|fe=o. 



(17) 



The quantity on the right hand side of this equation is 
plotted as a function of k 2 for I — 2 and I = 4 in Figure 
2. The k — > limit agrees well with the value obtained 
from the left hand side. 
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FIG. 2. Illustration of the Gubbins equation (17) for I — 2 
and 4 at the density p — 0.3o-g" 3 . Open symbols show right 
hand side of eqn. (17) vs. fc 2 , closed symbols the value at 
k — > obtained from the left hand side. The system size was 
N = 8000. 



III. RESULTS 

The orientational average of the pair correlation func- 
tions is shown in Figure 3 for different densities. It is 
basically given by their respective spherical harmonics 
coefficients with l,li,m,mi = 0. The pair distribution 
function exhibits several peaks at distances slightly 
larger than <7o, corresponding to the nearest neighbor 
shell, the next nearest neighbor shell etc. The peaks are 
less pronounced, but still present in the total correlation 
function. They disappear in the DCF. At low density, 
the DCF vanishes at distances beyond k<tq. At large 
densities, however, one observes broad oscillations which 
persist further. Oricntationally averaged, the structure 
of the DCF does not appear to be much simpler than 
that of the pair distribution function. 
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FIG. 3. Orientational average of the pair distribution func- 
tion p^ (a), the total correlation function h (b), and the 
direct correlation function c (c) vs. distance r for different 
densities. Data for different system sizes (p = 0.3cr~ 3 ) lie on 
top of each other. 

The qualitative difference between the DCF and the 
total correlation function becomes apparent when looking 
at coefficients which reflect the elasticity of the nematic 
fluid. Figure 4 shows as an example the coefficients with 
l\ = l 2 = I = 2, mi = — m 2 = 1 and m = 0. In the 
nematic phase, they decay at long distance like 1 /r, both 
for the pair distribution function and the total pair cor- 
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FIG. 4. Expansion coefficient with U = I = 2 and 
mi = — m,2 = 1 of the pair distribution function p' 2 ' (a), 
the total correlation function h (b), and the direct correla- 
tion function c (c) vs. distance r for the density p = 0.24o-^ 3 
(long dashed line) and p = 0.3cr^" 3 (solid line: system size 
N = 8000, dotted line: N = 4000, short dashed line: 
N = 1000). Dashed dotted line in the plot for h indicates 
extrapolation with a 1/r behavior. Inset shows data for h as 
a function of 1/r. 

relation function. The DCF no longer exhibits such a 
long ranged tail. Thus the DCF is indeed a quantity 
which characterizes the structure of nematic fluids on a 
local scale. 
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FIG. 5. Expansion coefficient with U = 2 and 
/ = m,i = m = of the pair distribution function p^ (a), 
the total correlation function h (b), and the direct correla- 
tion function c (c) vs. distance r for different densities. Dot 
dashed line in the plot for p^ indicates infinite range limit 
for p — 0.2S7ag 3 . Inset of (b) shows blowup of data for large 
distances. 

Next wc investigate the isotropic/nematic transition in 
more detail. Figure 5 illustrates how oricntational order 
emerges and gives rise to a constant long range contri- 
bution to the pair distribution function coefficients with 
li =/= 0, m, = I = m = 0. The constant contribution 
is subtracted in the total correlation function (cf. eqn. 
(6)). Close to the transition, however, the total correla- 



tion function features a slowly decaying, quasi long range 
tail (Fig. 5 b). It stems from near critical fluctuations 
with large correlation lengths close to the transition, due 
to the fact that the transition is only weakly first order. 
Like the elastic tails, these quasicritical long range tails 
are no longer present in the DCF. 

From the results shown so far (Figs 3, 4, and 5), the 
effect of nematic ordering on the DCF is not obvious. 
The influence of nematic symmetry breaking on correla- 
tion functions is best studied in a molecular-fixed frame 
system, which takes full advantage of the symmetries in 
the isotropic phase and thus allows to assess directly the 
loss of symmetries in the nematic phase. Following Gray 
and Gubbins 3 , we introduce a local coordinate system in 
which the z axis points in the direction 1*12 of the inter- 
molecular distance vector ( "molecular frame" , see figure 
6). The other two axes are chosen such that the x axis 
lies in the plane spanned by and the director n 45 . 
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FIG. 6. Illustration of the coordinate frames. Thick solid 
lines show coordinate axes of director frame, where the z axis 
points in the direction of the director n. Thick dashed lines 
indicate axes of the molecular frame, where the 2-axis points 
along the vector ri2 which connects the two particles. The x 
axis of the molecular frame lies in the plane which is spanned 
by ri2 and n. 

This frame is now used in the spherical harmonics ex- 
pansion (9) of functions which depend on particle orien- 
tations Ui. The intcrmolccular distance vector ri2 itself 
is still represented in the director frame. In an isotropic 
fluid, pair correlation functions must be rotationally in- 
variant, and all expansion coefficients except those with 
I = m = vanish. In a nematic fluid, this is no longer 
true 48 . Rotational invariance is broken if coefficients with 
I 7^ become nonzero. One would certainly expect this 
to happen for the pair distribution function and the total 
pair correlation function. In the case of the DCF, the sit- 
uation is less clear. Zhong and Petschek 32 ' 35 have argued 
that the DCF should remain rotationally invariant in an 
aligned fluid, since it reflects effective pair potentials be- 
tween particles. 
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FIG. 7. Expansion coefficient with £» = I = 2 and 
rrii = m = of the direct correlation function in the molec- 
ular frame vs. r for different densities. In the isotropic 
phase, this coefficient must vanish for symmetry reasons. At 
p — 0. 283(7,7 3 , a small nonzero remnants are observed due to 
finite size effects. 
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FIG. 8. Expansion coefficient with U — 2 and 
/ = rrii = m = of the direct correlation function in the 
molecular frame vs. r for different densities. 

An equation which converts expansion coefficients from 
the director frame into the molecular frame is derived in 
the appendix. Such conversion formulae can be found 
for isotropic fluids in the literature 3 . Our general case is 
more involved. 

Figure 7 shows an example of an expansion coefficient, 
C202020! which vanishes in the case of rotational invari- 
ance. In the isotropic phase, the values are indeed close 
to zero. In the nematic phase, a nonzero DCF coeffi- 
cient not only emerges, it may even grow quite large. At 
p = 0.3<Jq 3 the values are comparable to those of a sym- 
metry preserving coefficient of similar order (8). Hence 
the DCF reflects the broken symmetry in the nematic 



phase. Zhong's and Petschek's seemingly reasonable as- 
sumption is not valid for fluids of ellipsoids. 

In contrast, the symmetry preserving coefficient C202000 
seems barely affected by the transition (Fig. 8). The 
curves for p = 0.283(7^ 3 and p = 0.287<7(7 3 , slightly below 
and above the transition, lie almost on top of each other. 

It is instructive to compare the simulation results for 
the DCF with the results from the Percus-Yevick the- 
ory 4 . In the Percus-Yevick approximation, the direct cor- 
relation function is determined by the Ornstcin-Zernikc 
equation (7) and the additional closure relation 

c(ui,u 2 ,r 12 ) - (l-e v ^' kBT ) (/i(ui,u 2 ,r 12 ) + l), (18) 

where V\2 is the pair potential between the particles 1 and 
2. In a homogeneous isotropic fluid, this defines a closed 
set of equations. In an orientationally fluid, one needs 
an additional condition for the one particle distribution 
p^\u), e. g., eqn. (16). 

We have calculated the DCF in Percus-Yevick approx- 
imation for the density values of our simulation. We fol- 
lowed the procedure described by Letz and Latz 13 , except 
that we imposed eq. (16) instead of isotropy. Neverthe- 
less, we only obtained isotropic solutions. This is con- 
sistent with an argument of Zhong and Petschek, who 
showed that the diagrammatic expansion of the Ward 
identity, an equation related to (16), must contain di- 
agrams in the anisotropic case which are excluded by 
the Percus-Yevick closure 32 . Thus the Percus Yevick- 
closure is probably not compatible with spontaneous ori- 
cntational ordering. 

In the Percus-Yevick calculations, the spherical har- 
monics expansion was extended up to ^ max = 4. In addi- 
tion, we have also performed systematic calculations with 
^max = 2. Here we observed an interesting scenario: Over 
a wide range of densities, two solutions exist that meet in 
a bifurcation at an upper critical density p c and disappear 
beyond p c . One of the branches merges into the Onsager 
solution at low densities. The other yields much more 
structured correlation functions and is clearly unphysical. 
Unfortunately, we have not been able to analyze system- 
atically the situation at higher cutoff Z max . Our results 
suggest that the scenario at Z max = 4 might be similar. 
(For example, we found two solutions at p — 0.29cr^" 3 ). 
Therefore it may well be that the behavior observed at 
''max = 2 is a generic-feature of the Percus Yevick ap- 
proximation and not an artefact of the truncation of I. 
This would explain how the Percus-Yevick solution can 
disappear at high densities without exhibiting a nematic 
instability 18 . 

Figure 9 compares the Percus-Yevick prediction for 
the DCF coefficient with li — 2,1 = m = rrii = 
with the simulation data for densities p — 0.24er~ 3 and 
p = 0.283cr~ 3 . The agreement is good even for the higher 
density, which is just below the transition. Above the 
transition, we can contrast the DCF with the Percus- 
Yevick solution for the isotropic fluid of the same density. 
For example, the Percus-Yevick DCF for p = 0.3cr^" 3 , 
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shown in figure 9 b) , can be compared with the real DCF, 
shown in figure 8. The agreement is not very good. The 
real, nematic DCF exhibits much more structure than 
the Pcrcus-Yevick result. 




2,4 6 

r/G 

FIG. 9. Expansion coefficient with U = 2 and 
I = mi = m = of the direct correlation function in the 
molecular frame vs. distance r in Percus-Yevick approxima- 
tion and in simulations at the density p = 0.24(7,7 3 (a) and 
p = 0.283(7,7 3 (b). Simulation data are the same as in fig. 
8. (same data as figure 8). Dashed line shows Percus-Yevick 
prediction for p — 0.3(7,7 3 for comparison. 

IV. SUMMARY AND CONCLUSIONS 

To summarize, we have studied the local structure of 
fluids of uniaxial ellipsoids at different densities in the 
vicinity of the nematic/isotropic phase transition. In 
particular, we have calculated and analyzed the DCF. 
We find that the DCF suitably characterizes the local 
structure, in the sense that it is short ranged both in 
the isotropic and the nematic phase. This reflects the 
short range nature of the interactions in our model. In 
other respect, however, the DCF does not reproduce the 
properties of the pair interaction potential. Noteably, it 



loses rotational invariance in the nematic phase and fea- 
tures a substantial symmetry breaking contribution. In 
contrast, the symmetry preserving part of the DCF does 
not change dramatically at the nematic/isotropic transi- 
tion. 

Our results should be useful guides for future develop- 
ments of density functional theories for anisotropic flu- 
ids 49 , or of improved liquid state theories for nematic 
liquids in general. Moreover, we can use the DCF di- 
rectly to predict the structure of inhomogeneous fluids 
of ellipsoids within simple density functional or integral 
equation approaches 1 . Such work is currently under way. 
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APPENDIX 

In this appendix, we derive an equation which allows 
to convert expansion coefficients of pair correlation func- 
tions from the director frame to the molecular frame. 

We denote by F dir (ui, u 2 , r) and ^ mo '(ui, u' 2 , r) the 
representation of a correlation function F in the director 
frame and the molecular frame, respectively. These are 
expanded in spherical harmonics according to eqn. (9), 
hence 

F*> l5 u 2 ,r)= Y, F^m.imir) (19) 

l lt l 2 ,l 
m 1 ,m 2 ,m 

Y (imi (ui) Y hm2 (u 2 ) Y lm (r), 

and 

^(uijy^Kjy^K). (20) 

Now let lZ(fl) = 1Z(9, (j), x) be the rotation operator 
which carries the director frame into coincidence with 
the molecular frame (see figure 6), 

F mo ViX 2 ,r) = F* r (ftui,fti4r). (21) 

The angles 9 and 4> are the polar coordinates of the inter- 
molecular distance vector in the director frame, and 
X = 7r. The coefficients of F in the director frame and in 
the molecular frame are then related to each other by 
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pmol ( \ _ K l' 1 m' 1 l' 2 m' 2 V ' m' p dir ( \ ^[m'^m^l' m' _ ? r r 

r l' 1 m l 1 l' 2 m' 2 l'm'V ) — Iy l imi l 2 m 2 lm r l 1 m 1 l 2 m 2 lmV h "imnijmalra ~~ < J l 1 l 1 < J l 2 l 2 < J m 1 + 



m 2 +m,0 



(30) 



(22) 



KZhZ™ = J durduidi^^Cuijy^K) (23) 



x^ m ,(u r )H imi (Kui)y, 2ra2 (Ku' 2 )y lra (u r ). 
Next we apply the transformation formula 3 

H imi (^ui) - ^B^ffiJ'^MK), (24) 



M 



where D l ^ iM (fl) is the rotation matrix of the operator 
1Z, given by 



DL(^,^x)=e-^^(-l) fe (JyC i 



Irnnk 



(25) 



x (1 + a) ,-fc+I Tr 1 (l - a) fe " I!1 ^e- <nx 
with a = cos 9 and 



^ + m)!(i-m)!(/ + n)!(/-n)! 
° im " fc (J + m-A)!(i-n-fc)!fc!(fc-ro + n)!" [ ' 

The sum over fc is taken over values such that the argu- 
ments of the factorials in the denominator of Ci mn k are 
positive, i.e., max(0,m — n) < k < min(Z — n, I + m). 
Eq. 23 then becomes 



K2ZZ™ = / du r Y{! m ,{n r )Y lm (ur) (27) 



DZn^r^vDZn^rs,^. 

Using eqn. (25) with \ = tt and the relation 



1 (21 + 1) (I - m y 
Ylm -]j An (l + m)\ e PlM (28) 

(Pim is the associated Legendre polynomial), we obtain 

ri -l 1 mil 2 m 2 lm ~ hl' 1 hl' 2 m 1 +m 2 +m,m' \ / 

x(-) v(2^+i)(2z + i )y a+m)!(/ , +m0 , 

V \ f _1 ^fcl+fe/^ fl 

A / V Wimim'jti { -^l 2 m 2 m' 2 k 2 

ki ,k 2 
/■l 



1 



x - / daP llm ,{a)Pi m {a) (1 + a)' 1 



(— — m2+mi+m2+2fci+2/c2) 



x V(l-a)/(l + a) 



In a uniaxial nematic phase, molecular frame coefficients 
with m' 7^ must vanish. Thus eqn. (29) reduces to 



X ^ ' (~^) ^~ / 'l 1 m' 1 m 1 k 1 ^l 2 m' 2 m 2 k 2 

ki ,k 2 

i r 1 

x -y daiV(«) *W«) (l + a)' 1+i2 . 

x v /(l-«)/(l + «) 2fcl+2fe2+m ' 1+ " 4 (-l) m ' 1+ ^ 

In an isotropic phase, one can use this result with V = 
m' — and = — m' 2 = fh. 

The matrix for the inverse transformation from the 
molecular frame to the director frame the same, eqn. 
(30). We have performed several numerical tests to check 
this result: Pair distribution functions were determined 
from the simulation data both in the molecular frame and 
in the director frame. Then we applied the transforma- 
tion formula and checked that the result is the same 50 . In 
another test, we have transformed correlation functions 
into the molecular frame and back again, and compared 
the result with the original data. 
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